package com.li.fluids.sph
{
import com.li.fluids.sph.MpmMaterial;

public class MpmFluid
{
    private var _gsizeX:int;
    private var _gsizeY:int;
    private var _grid:Vector.<Vector.<MpmNode>>;
    private var _activeNodes:Vector.<MpmNode> = new Vector.<MpmNode>();
    private var _particles:Vector.<MpmParticle> = new Vector.<MpmParticle>();
    private var _points:Vector.<Number> = new Vector.<Number>();
    private var _cellWidth:Number;
    private var _cellHeight:Number;

    public function MpmFluid(width:Number, height:Number, cellWidth:Number, cellHeight:Number)
    {
        _cellWidth = cellWidth;
        _cellHeight = cellHeight;
        _gsizeX = width/cellWidth;
        _gsizeY = height/cellHeight;

        // Initialize grid.
        var i:int, j:int;
        _grid = new Vector.<Vector.<MpmNode>>(_gsizeX, true);
        for(i = 0; i < _gsizeX; i++)
        {
            _grid[i] = new Vector.<MpmNode>(_gsizeY, true);
            for(j = 0; j < _gsizeY; j++)
                _grid[i][j] = new MpmNode();
        }
    }

    public function get points():Vector.<Number>
    {
        return _points;
    }

    public function addParticleBlock(material:MpmMaterial, startx:Number, starty:Number, widthx:Number, widthy:Number):void
    {
        var i:uint, j:uint;
        for(i = startx; i < widthx; i++)
            for(j = starty; j < widthy; j++)
                addParticle(material, i + 4, j + 4, 0.0, 0.0);
    }

    public function addParticle(material:MpmMaterial, x:Number, y:Number, u:Number, v:Number):void
    {
        _particles.push(new MpmParticle(material, x, y, u, v));
        _points.push(_cellWidth*x, _cellHeight*y);
    }
    
    public function simulate():void
    {
        var n:MpmNode, p:MpmParticle;
        for each (n in _activeNodes)
            n.clear();
        _activeNodes.length = 0;

        var i:int, j:int;
        var x:Number, y:Number, phi:Number;
        var fx:Number = 0.0, fy:Number = 0.0;
        for each (p in _particles)
        {
            p.cx = int(p.x - 0.5);
            p.cy = int(p.y - 0.5);

            x = p.cx - p.x;
            p.px[0] = (0.5*x*x + 1.5*x + 1.125);
            p.gx[0] = (x + 1.5);
            x += 1.0;
            p.px[1] = (-x*x + 0.75);
            p.gx[1] = (-2.0*x);
            x += 1.0;
            p.px[2] = (0.5*x*x - 1.5*x + 1.125);
            p.gx[2] = (x - 1.5);

            y = p.cy - p.y;
            p.py[0] = (0.5*y*y + 1.5*y + 1.125);
            p.gy[0] = (y + 1.5);
            y += 1.0;
            p.py[1] = (-y*y + 0.75);
            p.gy[1] = (-2.0*y);
            y += 1.0;
            p.py[2] = (0.5*y*y - 1.5*y + 1.125);
            p.gy[2] = (y - 1.5);

            for(i = 0; i < 3; i++)
            {
                for(j = 0; j < 3; j++)
                {
                    n = _grid[p.cx + i][p.cy + j];
                    if(!n.active)
                    {
                        _activeNodes.push(n);
                        n.active = true;
                    }
                    phi = p.px[i]*p.py[j];
                    n.m += phi*p.mat.m;
                    n.d += phi;
                    n.gx += p.gx[i]*p.py[j];
                    n.gy += p.px[i]*p.gy[j];
                }
            }
        }

        var density:Number, pressure:Number, weight:Number;
        var n01:MpmNode, n02:MpmNode;
        var n11:MpmNode, n12:MpmNode;
        var cx:int, cy:int;
        var cxi:int, cyi:int;

        var pdx:Number, pdy:Number;
        var C20:Number, C02:Number, C30:Number, C03:Number;
        var csum1:Number, csum2:Number;
        var C21:Number, C31:Number, C12:Number, C13:Number, C11:Number;

        var u:Number, u2:Number, u3:Number;
        var v:Number, v2:Number, v3:Number;

        for each (p in _particles)
        {
            cx = int(p.x);
            cy = int(p.y);
            cxi = cx + 1;
            cyi = cy + 1;

            n01 = _grid[cx][cy];
            n02 = _grid[cx][cyi];
            n11 = _grid[cxi][cy];
            n12 = _grid[cxi][cyi];

            pdx = n11.d - n01.d;
            pdy = n02.d - n01.d;
            C20 = 3.0*pdx - n11.gx - 2.0*n01.gx;
            C02 = 3.0*pdy - n02.gy - 2.0*n01.gy;
            C30 = -2.0*pdx + n11.gx + n01.gx;
            C03 = -2.0*pdy + n02.gy + n01.gy;
            csum1 = n01.d + n01.gy + C02 + C03;
            csum2 = n01.d + n01.gx + C20 + C30;
            C21 = 3.0*n12.d - 2.0*n02.gx - n12.gx - 3.0*csum1 - C20;
            C31 = -2.0*n12.d + n02.gx + n12.gx + 2.0*csum1 - C30;
            C12 = 3.0*n12.d - 2.0*n11.gy - n12.gy - 3.0*csum2 - C02;
            C13 = -2.0*n12.d + n11.gy + n12.gy + 2.0*csum2 - C03;
            C11 = n02.gx - C13 - C12 - n01.gx;

            u = p.x - cx;
            u2 = u*u;
            u3 = u*u2;
            v = p.y - cy;
            v2 = v*v;
            v3 = v*v2;
            density = n01.d + n01.gx*u + n01.gy*v + C20*u2 + C02*v2 + C30*u3 + C03*v3 + C21*u2*v + C31*u3*v + C12*u*v2 + C13*u*v3 + C11*u*v;

            pressure = density - 1.0;
            if(pressure > 2.0)
                pressure = 2.0;

            fx = 0.0;
            fy = 0.0;

            if(p.x < 4.0)
                fx += p.mat.m*(4.0 - p.x);
            else if(p.x > _gsizeX - 5)
                fx += p.mat.m*(_gsizeX - 5 - p.x);

            if(p.y < 4.0)
                fy += p.mat.m*(4.0 - p.y);
            else if(p.y > _gsizeY - 5)
                fy += p.mat.m*(_gsizeY - 5 - p.y);

            for(i = 0; i < 3; i++)
            {
                for(j = 0; j < 3; j++)
                {
                    n = _grid[(p.cx + i)][(p.cy + j)];
                    phi = p.px[i]*p.py[j];
                    n.ax += -((p.gx[i]*p.py[j])*pressure) + fx*phi;
                    n.ay += -((p.px[i]*p.gy[j])*pressure) + fy*phi;
                }
            }
        }

        for each (n in _activeNodes)
        {
            if(n.m > 0.0)
            {
                n.ax /= n.m;
                n.ay /= n.m;
                n.ay += 0.03;
            }
        }

        var mu:Number, mv:Number;
        for each (p in _particles)
        {
            for(i = 0; i < 3; i++)
            {
                for(j = 0; j < 3; j++)
                {
                    n = _grid[(p.cx + i)][(p.cy + j)];
                    phi = p.px[i]*p.py[j];
                    p.u += phi*n.ax;
                    p.v += phi*n.ay;
                }
            }
            mu = p.mat.m*p.u;
            mv = p.mat.m*p.v;
            for(i = 0; i < 3; i++)
            {
                for(j = 0; j < 3; j++)
                {
                    n = _grid[(p.cx + i)][(p.cy + j)];
                    phi = p.px[i]*p.py[j];
                    n.u += phi*mu;
                    n.v += phi*mv;
                }
            }
        }

        for each (n in _activeNodes)
        {
            if(n.m > 0.0)
            {
                n.u /= n.m;
                n.v /= n.m;
            }
        }

        var gu:Number, gv:Number;
        var particleIndex:uint;
        for each (p in _particles)
        {
            gu = 0.0;
            gv = 0.0;
            for(i = 0; i < 3; i++)
            {
                for(j = 0; j < 3; j++)
                {
                    n = _grid[(p.cx + i)][(p.cy + j)];
                    phi = p.px[i]*p.py[j];
                    gu += phi*n.u;
                    gv += phi*n.v;
                }
            }
            p.x += gu;
            p.y += gv;
            p.u += 1.0*(gu - p.u);
            p.v += 1.0*(gv - p.v);
            if(p.x < 1.0)
            {
                p.x = (1.0 + Math.random()*0.01);
                p.u = 0.0;
            }
            else if(p.x > _gsizeX - 2)
            {
                p.x = (_gsizeX - 2 - Math.random()*0.01);
                p.u = 0.0;
            }
            if(p.y < 1.0)
            {
                p.y = (1.0 + Math.random()*0.01);
                p.v = 0.0;
            }
            else if(p.y > _gsizeY - 2)
            {
                p.y = (_gsizeY - 2 - Math.random()*0.01);
                p.v = 0.0;
            }

            // Update output positions.
            var pointsIndex:uint = 2*particleIndex;
            _points[pointsIndex] = p.x*_cellWidth;
            _points[pointsIndex + 1] = p.y*_cellHeight;
            particleIndex++;
        }
    }
}
}
